library(readxl)
library(ggplot2)
library(dplyr)
library(scales)
library(tidyr)
library(here)
library(utils)
library(sandwich)
library(lmtest)
library(easystats)
library(ggeffects)
library(vtable)
library(fixest)
library(MASS)

rm(list=ls(all = TRUE))

#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
# Load data
#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
main1_q <- readRDS(here("main1_q.rds"))
main2_q <- readRDS(here("main2_q.rds"))
main1_b <- readRDS(here("main1_b.rds"))
main2_b <- readRDS(here("main2_b.rds"))

table3_q <- readRDS(here("table3_q.rds"))
table3_b <- readRDS(here("table3_b.rds"))

parallel <- readRDS(here("parallel.rds"))
newcomer <- readRDS(here("newcomer.rds"))
questions <- readRDS(here("questions.rds"))
district <- readRDS(here("district.rds"))
reform <- readRDS(here("reform.RDS"))


#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
# Main text
#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
#########
# TABLE 1
#########
table1 <- main1_q %>% 
  filter(post == 0)

table1$age <- as.numeric(table1$age)
table1$leg_office <- as.numeric(table1$leg_office)
table1$localexp2 <- as.numeric(table1$localexp2)

# summary statistics with ttest
st(table1,
   group = 'maj_tier',
   vars = c('age', 'male', 'gov_office', 'leg_office', 'local', 'localexp2'),
   out = 'viewer',
   group.test = TRUE,
   digits = 2,
   fixed.digits = T)

#########
# TABLE 2
#########
table2_q <- main1_q %>% # use main2_q for post == 2
  filter(post == 0) # change post == 0/1/2 to get different rows
st(table2_q,
   group = 'maj_tier',
   vars = c('qreg_share'),
   out = 'viewer',
   group.test = F,
   digits = 2,
   fixed.digits = T)

table2_b <- main1_b %>% # use main2_b for post == 2
  filter(post == 0) # change post == 0/1/2 to get different rows
st(table2_b,
   group = 'maj_tier',
   vars = c('breg_share'),
   out = 'viewer',
   group.test = F,
   digits = 2,
   fixed.digits = T)

#########
# TABLE 3 
#########
# Questions
### SMD MPs
SMD_01_q <- table3_q %>% 
  filter(maj_tier == 1)
sumtable(SMD_01_q, 
           group = 'reelected', 
           vars = 'qreg_share',
           out = 'viewer',
           group.test = T, 
           digits = 2,
           fixed.digits = T)

### MMD MPs
MMD_01_q <- table3_q %>% 
  filter(maj_tier == 0)
sumtable(MMD_01_q, 
         group = 'reelected', 
         vars = 'qreg_share',
         out = 'viewer',
         group.test = T, 
         digits = 2,
         fixed.digits = T)

# Bills
### SMD MPs
SMD_01_b <- table3_b %>% 
  filter(maj_tier == 1)
sumtable(SMD_01_b, 
         group = 'reelected', 
         vars = 'breg_share',
         out = 'viewer',
         group.test = T, 
         digits = 2,
         fixed.digits = T)

### MMD MPs
MMD_01_b <- table3_b %>% 
  filter(maj_tier == 0)
sumtable(MMD_01_b, 
         group = 'reelected', 
         vars = 'breg_share',
         out = 'viewer',
         group.test = T, 
         digits = 2,
         fixed.digits = T)

#########
# TABLE 4 
#########
### (1)
reg1_q <- lm(qreg_share ~ post*as.factor(maj_tier), data = main1_q)
coeftest(reg1_q, vcov = vcovCL, cluster = ~ persona)

### (2)
reg2_q <- lm(qreg_share ~ post*as.factor(maj_tier) +  as.factor(gov_office) + as.factor(leg_office) + as.factor(local)
                 + party_merge, 
                 data = main1_q)
coeftest(reg2_q, vcov = vcovCL, cluster = ~ persona)

### (3)
reg3_q <- lm(qreg_share ~ as.factor(post)*as.factor(maj_tier), data = main2_q)
coeftest(reg3_q, vcov = vcovCL, cluster = ~ persona)

### (4)
reg4_q <- lm(qreg_share ~ as.factor(post)*as.factor(maj_tier) +  as.factor(gov_office) + as.factor(leg_office) + as.factor(local)
             + party_merge, 
           data = main2_q)
coeftest(reg4_q, vcov = vcovCL, cluster = ~ persona)

##########
# FIGURE 4 
########### 
main2_q$maj_tier <- as.character(main2_q$maj_tier)
main2_q$post <- as.character(main2_q$post)
main2_q$leg_office <- as.numeric(main2_q$leg_office)
reg2_q <- lm(qreg_share ~ post*as.factor(maj_tier) +  as.factor(gov_office) + as.factor(leg_office) + as.factor(local)
             + party_merge, 
             data = main2_q)

fig4 <- ggpredict(reg2_q, terms = c("post","maj_tier"), ci.lvl = 0.9, 
                   condition = c(leg_office = "1", local = "1", gov_office ="0"))

# To make results comparable, only MPs re-elected in the third term were used 
# for Figure 4a. 
plot(fig4, connect.lines = F, colors = 'social') + 
  labs(
    x = "Period", 
    y = "% of regional questions", 
    title = "") +
  labs(title = NULL, colour = "Tier")

#########
# TABLE 5 
#########
### (1)
reg1_b <- lm(breg_share ~ post*as.factor(maj_tier), data = main1_b)
coeftest(reg1_b, vcov = vcovCL, cluster = ~ persona)

### (2)
reg2_b <- lm(breg_share ~ post*as.factor(maj_tier) +  as.factor(gov_office) + as.factor(leg_office) + as.factor(local)
             + party_merge, 
             data = main1_b)
coeftest(reg2_b, vcov = vcovCL, cluster = ~ persona)

### (3)
reg3_b <- lm(breg_share ~ as.factor(post)*as.factor(maj_tier), data = main2_b)
coeftest(reg3_b, vcov = vcovCL, cluster = ~ persona)

### (4)
reg4_b <- lm(breg_share ~ as.factor(post)*as.factor(maj_tier) +  as.factor(gov_office) + as.factor(leg_office) + as.factor(local)
             + party_merge, 
             data = main2_b)
coeftest(reg4_b, vcov = vcovCL, cluster = ~ persona)

#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
# Appendix
#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

###########
# FIGURE C1 
###########
main2_q <- readRDS("main2_q.rds")
reg2_q <- lm(qreg_share ~ post*as.factor(maj_tier) +  as.factor(gov_office) + as.factor(leg_office) + as.factor(local)
             + party_merge, 
             data = main1_q)

eqtest <- equivalence_test(reg2_q, range = c(-0.12, 0.12), ci = 0.90, rule = "classic") # change to [-0.13, 0,13] to see that HO becomes accepted
plot(eqtest)

###########
# FIGURE D1 
###########
ggplot(parallel, aes(x = term, y= geo_rep, group=maj_tier)) +
  geom_line(aes(color = maj_tier))+
  geom_point(aes(color= maj_tier)) +
  ylab("Geographical representation") +
  xlab("") +
  theme_minimal()


###########
# FIGURE E1 
###########
reg_newcomer <- lm(qreg_share ~ experienced*maj_tier +  gov_office + leg_office + local + party_1, 
                   data = newcomer)
plot_newcomer <- ggpredict(reg_newcomer, terms = c("experienced","maj_tier"), ci.lvl = 0.9)

plot(plot_newcomer, connect.lines = F, colors = 'social') + 
  labs(
    x = "Experience", 
    y = "% of regional questions", 
    title = "") +
  scale_x_continuous(labels = c("New entrant", "MP with prior experience"), breaks = c(0, 1)) +
  labs(title = NULL, colour = "Tier")


###########
# TABLE B1 
###########
logit1 <- glm(reg_qs ~ post*maj_tier, 
                family = binomial(link = "logit"), 
                data = questions)
coeftest(logit1, vvcov = vcovCL, cluster = ~ persona)


logit2 <- glm(reg_qs ~ post*maj_tier +  gov_office + leg_office + local 
                + party_merge, 
                family = binomial(link = "logit"), 
                data = questions)
coeftest(logit2, vvcov = vcovCL, cluster = ~ persona)


###########
# TABLE B2 
###########
### FIxed effects models
main1_q <- main1_q %>%
  mutate(treated = maj_tier == 1 & post == 1)

fixed1 <- feols(qreg_share ~ treated | maj_tier + post, data = main1_q)
summary(fixed1, cluster = ~ persona)

fixed2 <- feols(qreg_share ~ treated + gov_office + leg_office + local 
                + party_merge | maj_tier + post  , data = main1_q)
summary(fixed2)

fixed3 <- feols(qreg_share ~ treated + gov_office + leg_office + local 
                + party_merge | maj_tier + post + persona  , data = main1_q)
summary(fixed3, cluster = ~ persona)


###########
# TABLE B3 
###########
regref1 <- lm(qreg_share ~ post*as.factor(maj_tier), data = reform)
coeftest(regref1, vcov = vcovCL, cluster = ~ persona)

regref2 <- lm(qreg_share ~ post*as.factor(maj_tier) +  as.factor(gov_office) + as.factor(leg_office) + as.factor(local)
             + party_merge, 
             data = reform)
coeftest(regref2, vcov = vcovCL, cluster = ~ persona)


###########
# TABLE B4 
###########
### Count models
count1 <- glm.nb(count_qreg ~ post*maj_tier + offset(log(days_office)), 
                 control=glm.control(maxit=50), data = main1_q)
coeftest(count1, vcov = vcovCL(count1, cluster = ~ persona, type = "HC1"))


count2 <- glm.nb(count_qreg ~ post*maj_tier + gov_office + leg_office + local 
                  + party_merge + offset(log(days_office)),
                  control=glm.control(maxit=50),
                  data = main1_q)
coeftest(count2, vcov = vcovCL(count2, cluster = ~ persona, type = "HC1"))



###########
# TABLE B5 
###########
### Count controlling for total questions
counttot1 <- glm.nb(count_qreg ~ post*maj_tier + count_qtot + offset(log(days_office)), 
                 control=glm.control(maxit=50), data = main1_q)
coeftest(counttot1, vcov = vcovCL(counttot1, cluster = ~ persona, type = "HC1"))


counttot2 <- glm.nb(count_qreg ~ post*maj_tier + gov_office + leg_office + local 
                 + party_merge + count_qtot + offset(log(days_office)),
                 control=glm.control(maxit=50),
                 data = main1_q)
coeftest(counttot2, vcov = vcovCL(counttot2, cluster = ~ persona, type = "HC1"))


###########
# TABLE B6 
###########
nozero <- main1_q %>% 
  filter(count_qtot > 0)

regnz1 <- lm(qreg_share ~ post*as.factor(maj_tier), data = nozero)
coeftest(regnz1, vcov = vcovCL, cluster = ~ persona)

regnz2 <- lm(qreg_share ~ post*as.factor(maj_tier) +  as.factor(gov_office) + as.factor(leg_office) + as.factor(local)
             + party_merge, 
             data = nozero)
coeftest(regnz2, vcov = vcovCL, cluster = ~ persona)


###########
# TABLE B7 
###########
regdistr1 <- lm(qreg_share ~ post*maj_tier, 
                 data = district)
coeftest(regdistr1, vcov = vcovCL, cluster = ~ persona)


regdistr2 <- lm(qreg_share ~ post*maj_tier +  gov_office + leg_office + local 
                 + party_merge, 
                 data = district)
coeftest(regdistr2, vcov = vcovCL, cluster = ~ persona)

